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ABSTRACT 



We construct self-consistent dynamical models for disk galaxies with triax- 
ial, cuspy halos. We begin with an equilibrium, axisymmetric, disk-bulge-halo 
system and apply an artificial acceleration to the halo particles. By design, this 
acceleration conserves energy and thereby preserving the system's differential 
energy distribution even as its phase space distribution function is altered. The 
halo becomes triaxial but its spherically-averaged density profile remains largely 
unchanged. The final system is in equilibrium, to a very good approximation, 
so long as the halo's shape changes adiabatically. The disk and bulge are "live" 
while the halo is being deformed; they respond to the changing gravitational 
potential but also influence the deformation of the halo. We test the hypothesis 
that halo triaxiality can explain the rotation curves of low surface brightness 
galaxies by modelling the galaxy F568-3. 

Subject headings: Galaxy: kinematics and dynamics — methods: statistical - 
methods: N-body simulations — cosmology: dark matter 



INTRODUCTION 



Dark matter halos - at least the ones found in cosmological simulations - have a num- 
ber of universal traits. Most famously, their density profiles hav e a shape that is nearly 
indep endent of mass, formation epoch, and cos mological model (INavarro. Frenk fc White 
996 ). Their angular m omentum distribution (IBullock et al.l 12001 ) , pha se space density 



(ITaylor fc Navarro! 1200 ll ). and velocity anisotropy ([Hansen fc Moord 120061 ) profiles also ap 



pear to follow universal forms. In other respects, halos are rather diverse. Simulated halos 
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are ty p ically triaxial with a xis ratios that rang e from 0.6 to 1 (Se e iDubinski fc Carlberg 
( 1l99lh ; I Warren et al.l (119921 ) and more recently, iNovak et al.l (120061 )). The shapes of real 
halos are more difficult to determine but promising obser vational approaches do exist. 
Probes of the Galactic halo include flaring of the gas di sk (jOlling fc Merrifieldl |2000| ) and 
tidal streams of satellite galaxies (IJohnston et al.l Il999l ). The shapes of halos in other 
galaxies can be determined, at least statistically, by weak gravitational l ensing surveys 
feoekstra. Yee. fc Gladder^booi [Mandelbaum et alil2006l : lparker et alil2007h . On the other 
hand, triaxiality can bias attempts to determine a halo's den sity profile from the rotation 
curve of the disk that sits within it (lHayashi &: Navarro! 120061 ). 



Our main goal in writing this paper is to introduce a novel scheme to generate self- 
consistent dynamical models for disk galaxies with triaxial halos. Our models can be tailored 
to fit observational data for specific galaxies and therefore provide a testing ground to study 
the disk-halo connection. We consider the effects of halo triaxiality on the rotation curves of 
low surface brightness galaxies (LSBs) and briefly discuss other applications of the models. 

Halos in simulati ons of a cold dark matter (CD M) universe have central cusps with 
p oc r~ 7 where 7 ~ 1 (jNavarro. Frenk fc White! 119961 ) . In dark matter-dominated galaxies, 
this density profile would seem to imply a rotation curve where v oc r 1 / 2 as r — ► 0. By 
contrast, a halo with a constant density core implies v cx r as r — > 0. LSBs, which are 
believed to be dark-matter dominated at small radii, have rotation curves that generally favor 
a constant density core over a 7 = — 1 cusp. This result represents one of the most serious 
challe nges to the CDM scenario (lMoorelll994j ; iFlores fc Primacklll994 iMcGaugh k, de Blok 
19981 ) and has inspired some rather exotic alternatives. De Blok & McGaugh (1998), for 
example, sug gested th a t LSB rot ation curves could be explained by Modified Newtonian 
Gravity while iFirmanil (120001 ) and lMo fc Mad (120001 )) invoked dark matter self-interactions 
to flatten the central cusp of the halo. 

The connection described above between a galaxy's rotation curve and the intrinsic den- 
sity profile of its halo assumes that the halo is spherically symmetric. However, if a galaxy's 
halo is triaxial, then gas in the disk will move on non-circular orbits and, under certain con- 
ditions, the observed rotation c urve will rise approximate ly lin early even if the intrinsic halo 
densit y profile has a steep cusp. lHayashi &: Navarro! (120061 ) and lHayashi. Navarro, fc Springe! 
( 120071 ) presented this argument as a means of reconciling LSB rotation curves with the pre- 
dictions of the CDM model. 



Hayashi fc Navarro! (120061 ) and lHayashi. Navarro, fc Springell (120071 ) derived model ro 



tation curves by calculating closed orbits in the potential generated by a triaxial halo. In this 
paper, we derive rotation curves by making pseudo-observations of a disk that is embedded 
in the halo. Deviations from axial symmetry in disk and halo are generated concurrently 
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and self-consistently. 



A number of methods exist for construct ing; models of, and embedding disks in, tri- 
axial halos. For example, iMoore et al.l (120041 ) sho w that the r e mnan t of a major merger 
between two equilibrium spherical halos is triaxial. iBailin et al.l (120071 ) describe how to set 
up an equilibrium disk in a combined halo-disk potential. Our approach produces an N- 



body galaxy comp 
be included, as in 



ete with disk, bulge , and triaxial halo. (Central black holes may also 



Widrow fc Dubinskil (120051 ) .) It is inspired by the method outlined in 



Holley-Bockelmann et al.l (120011 ). In that scheme, dubbed "adiabatic squeezing", particles 
of an equilibrium halo are subjected to an artificial drag by modifying the force law of a 
standard N-body code and evolving the system forward in time. A triaxial halo is created 
if the drag has a different strength along three orthogonal directions, that is, along what 
become the three principle axes of the halo. The final model will be in equilibrium, to a good 
approximation, so long as the timescale for the halo's shape to change is slow as compared 
to the typical orbital timescale of the system. 

Adiabatic squeezing causes a halo to shrink in size. For an isolated halo, this shrink- 
ing can be reversed by simply rescaling the positions and velocities of the particles. Ob- 
viously, this method is unsuitable for disk-bulge-halo systems since the disk and bulge 
would be disrupted in an unphysical way. We propose a modification of this method in 
which drag is applied along one axis and "negative drag" is applied along the other two or 
vice versa depending on whether one wants a prolate or oblate halo. We require that for 
each particle, the change in energy due to the artificial drag force is zero. In this way, we 
change t he phase space distr i butio n of the particles but not their energy distribution. As 
noted in lBinney fc Tremaind (119871 ). if two systems have the same energy distribution, their 
spherically-averaged density profiles will be very similar even if their phase space distribution 
functions are different. 



Our starting point is the equilibrium model of IWidrow. Pym. &: Dubinskil (120071 ) which 
comprises a Sersic bulge, cuspy dark halo, and exponential disk. The model is described 
in terms of a phase space distribution function (DF) which, in turn, is a function of the 
integrals of motion. In the current version of the model, the halo component of the DF 
depends only on the energy. In the absence of a disk, the halo is spherically symmetric. 
With the disk included, the halo is flattened slightly but is still axisymmetric. Adiabatic 
deformation allows us to extend our disk-bulge-halo model to systems with triaxial halos. 

In Section 2, we describe our method and construct an example of an isolated triaxial 
halo. We then consider the LSB galaxy F568-3. In Section 3, we present axisymmetric, 
equilibrium models for this galaxy based on its published surface brightness profile and 
circular speed curve. In Section 4, we show how transforming the axisymmetric halo in one 
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of these models into a triaxial halo changes the shape of the rotation curve. We conclude 
in Section 5, by summarizing our results and briefly discussing further applications of the 
method. 



2. METHOD 

We begin with an N-body equilibrium halo and evolve the system forward in time 
using standard techniques augmented by an artificial, energy-conserving acceleration. To be 
precise, we introduce an acceleration into the equations of motion given by 

(ft ~ ft) vj + M m 

a x = f- v x , (1) 

v 

a y = ~2 v v ' ( 2 ) 

and 

a z = 5 ~v z (3) 

v 

where v is the speed of the pa rticle and a • v = 0, as required. The coefficients ft and ft 



are time-dependent. Following iHolley-Bockelmann et al.l (120011 ) we assume that ft and ft 



"turn on" at t = and increase to their respective maximum values over a period Tq with 
a time-dependence given by 

ft = ft, m a,(3(t/T G ) 2 -2(t/T G ) 3 ) . (4) 

ft and ft remain constant for a time Tq before decreasing to zero over a time Tjj. 

A few comments regarding the parameters ft, ft, T G , Tc, and To are in order. First, 
the degree by which the halo departs from spherical symmetry is given, roughly, by the 
integral J flidt = ft (Tq + T c + T D ). Rescaling the fts by a factor / and the T's by a factor 
f^ 1 leaves the final halo shape unchanged with one important caveat. sets the timescale 
over which the halo's shape changes and therefore must be longer than its dynamical time 
~ cih/<Jh in order to maintain adiabaticity. 

In the case of an isolate halo equations [Ul3] admit several discrete symmetries. For 
example, interchanging ft and ft is equivalent to interchanging x and y. Similar symmetries 
are listed in Table 1. The presence of a disk breaks these symmetries. In general, increasing 
ft causes the system to expand along the x-axis while increasing @ 2 causes the system to 
expand along the y-axis. Models with (3\ = ft are axisymmetric about the z-axis. 
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As an illustration, we transform an isolated, spherically-symmetric halo into one that is 
triaxial. We begin with a halo whose density profile is given by 

2 2 ' 7 ^ 1 . ( r-r h \ 

Phaio = -j — 2~ . . rz- — 3— erfc — =— . 5 

(r/a h y (1 + r/a h y 7 \V25r h J 

For this example, we set 7 = 1 (the NFW value), = fOkpc, o h = 100 km s -1 , r h = 100 kpc 
and 5r h = 10 kpc. We choose (3x, max = 0.24 Gyr" 1 , (3 2 ,max = 0.08 Gyr" 1 , T G = T D = 1 Gyr 
and t c = 3 Gyr. As required, the characteristic timescale for the halo, ah/(Th = 100 Myr, is 
short compared to the timescale, = 4 — 12 Gyr, associated with the artificial force. Our 
choice for the yields a halo whose short axis is along the ^-direction and whose long axis 
is along the x-direction. 



The model is evolved forward i n time using th e N-body code from IStifil (120031 ) which is 



based on the algorithm described in lDehnenl (12000I ). The code uses a multipole expansion for 



cell-cell interactions; computational costs scale approximately linearly with particle number 
N. The softening length is 200 pc and the timestep is 1 Myr. The system is evolved for a 
period of 15 Gyr. 

In Figure [TJ we show a contour plot of the projected surface density of the halo along 
the three principle axes. Note that the departure from spherical symmetry is strongest in 
the inner parts of the halo. In Figure [2j we show the axis ratios as a function of time. To 
be precise, we model the density field as an ellipsoidal distribution, 

y 2 ^ 

p = p(f) where f = x 2 + — H — - . (6) 

b l c 



The p arameters b and c are calculated through an iterative procedure as outlined in lDubinski fc Carlberg 



( 1l99ll ). We show the results using the inner third of the particles, the inner two thirds of the 
particles, and all of the particles. Again, we see that the halo is more spherical in the outer 
parts. Note that the axis ratios at all radii oscillate a bit at t — 5 Gyr, the time when the 
artificial acceleration is turned off. After this time, the axis ratios in the inner two thirds 
of the halo quickly settle down to constant values. The oscillations damp more slowly in 
the outer parts of the halo where the dynamical time is not much shorter than (3^ . One 
can minimize the oscillations by increasing the T's and decreasing the /3's but at the cost of 
additional computation time. 

In Figure [31 we show the spherically-averaged differential mass profile, dMj dr oc r 2 p 
for the initial model and for the deformed model at t = 6 Gyr and t = 12 Gyr and compare 
with equation 5. Also shown is the density profile calculated from the initial conditions 
and the density profile for the system evolved to 12 Gyr with no artificial acceleration. The 
former illustrates the role mass resolution plays on the measured density profile while the 
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latter illustrates the effects of force softening and two-body relaxation. We see that the 
spherically-averaged density profile is preserved to within the fluctuations introduced by 
these other effects. 

Table 2 presents results for the axes ratios for other choices of f3\ and /?2- Note that mod- 
els la-c (and likewise models 2a-b and models 3a-c) are equivalent through the symmetries 
described in Table 1. 



3. THE LSB GALAXY F568-3 

In this section, we construct axisymmetric, equilibrium models for F568-3, an LSB 
galaxy which has appeared in a number of studies. We describe our general axisymmetric 
disk-bulge-halo models, review published photometric and kinematic observations for this 
galaxy, and discuss the statistical techniques used to tailor the model to the data. 



3.1. Equilibrium Models for Disk- Bulge- Halo Systems 



Our starting point is the dynamical galactic model described in lWidrow. Pym. fc Dubinski 



( 120071 ). The model is axisymmetric and comprises an exponential disk, a Sersic bulge, and a 
halo whose density profile is given by equation [5l DFs for the bulge and halo are functions of 
the energy, E, and constructed via an Abel integral transform. The D F for the disk is con- 
struct ed from three integrals of motion following the method outlined in lKuijken fc Dubinski 
( 119951 ). The total DF for the composite system self-consistently satisfies the collisionless 
Boltzmann and Poisson equations. 

The bulge has a spherically-averaged density profile given, to a good approximation, by 



PbulgelH = Pb 



R r 



This density profile yields the Sersic law, 



Sfr) 



S e 



-fe(i?/J?e) 1/r ' 



(7) 



(8) 



for the projected mass density provided one setsp = 1— 0.6097/n+0.05563/n 2 ( iPrugniel fc Simien 
19971 ; iTerzic fc Graham! 120051 ). E , R e and n are free parameters while the constant b is ad- 
justed so that R e encloses half the total projected light or mass. 



The disk DF depends on E, the angular momentum about the symmetry axis, L ZJ 
and an approximate integral of motion, E z , which corresponds to the energy associated 
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with vertical motions of stars in the disk. The DF is adjusted so that the intrinsic three- 
dimensional density distribution and velocity dispersion profile are given, respectively, by 

p disk (R, z) = p exp-*'** sech 2 (z/z d )etfc ((R - R ont ) /6R ont ) . (9) 

and 

o 2 B {R)=a 2 m ^ V {-R/R a ) . (10) 



3.2. Surface Brightness Profile and Rotation Curve for F568-3 



Multi-band photometry for the LSB galaxy F568-3 is presented in lde Blok. van der Hulst. &: Bothun 

(119951 ) . The galaxy resembles a normal late-type galaxy exhibiting a disk and faint spiral 



arms. However, its c entral B-band s urface brightness is more than a magnitude fainter than 
the Freeman value ([Freeman! 1 19701 ) placing it squarely in the category of LSBs. For the 
pu rpose of modelling the galaxy, we use the R-band surface brightness profile from Figure 2 
of Jde Blok. van der Hulst. k Bothuij[l99~5} ). 



High-resolution rotation curves for F568-3 are described in lMcGaugh. Rubin, fc de Blok 



(120011 ) . The circular speed rises approximately linearly to 80 km s 1 within 4 kpc. It continues 
to rise beyond this radius reaching a maximum value of ~ 100 km s -1 at R ~ 12 kpc. 



3.3. Markov Chain Monte Carlo Analysis of F568-3 

We use Bayesian statistics and the Markov Chain Monte Carlo (MCMC) method to find 
suitable axisymmetric models for F568-3. MCMC provides an efficient means of mapping out 
the likelihood function over the full multi-dimensional parameter space and has a number of 
advantages over traditional maximum likelihood techniques. 

For a particular choice of model parameters, one can construct a likelihood function 
which quantifies the agreement between the model and the data. Maximization techniques, 
such as the simplex algorithm, allow one to hone in on the "best-fit" model. However, with 
a large number of parameters, the likelihood function may become difficult to characterize 
with many false maxima. Moreover, the computational costs of simple algorithms, such as 
grid-based searches, become prohibitive. 

The goal of our MCMC analysis is to calculate the posterior probability density function, 
p(M\D, I), of a Galactic model, M, given data, D, and prior information, /. From Bayes' 
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theorem we have 



p{M\D,I) 



p(M\I)p(D\M,I) 



(11) 



p(D\I) 

where p(M\I) is the prior probability density and p(D\I) = J dAp(M\D, I) is a normal- 
ization factor. In MCMC, one constructs a sequence or chain of models through parameter 
space chosen according to a prescribed algorithm. The distribution of models along the chain 
will be proportional to p (M\D, I) provided the chain is sufficiently long. 



In this work, we use the Metrop olis-Hastings algorithm ([Metropolis et al.lll953l ; Hastings 



1970l ) as outlined in iGregoryl (120051 ). The first model in the chain is chosen at random. A 
candidate for the second model is chosen by taking a step in parameter space according to a 
proposal distribution. Let 1Z be the ratio of the likelihood function of the candidate to that of 
the first model. The candidate is accepted a fraction, /, of the time where / = min{l,7?.}. 
Otherwise, the second model is taken to be identical to the first model. The process is 
repeated to find the third model and so forth. 

Care must be taken in selecting a proposal distribution. If the step size is too short, 
the chain moves slowly through parameter space and the time required to fully explore 
parameter space becomes prohibitively large. On the other hand, if the typical step size is 
too large, the acceptance rate wil l be very low. We use an iterative approach, as outlined in 
Widrow. Pym. fc Dubinskil (120071 ). to choose an efficient proposal distribution. 



Since the data used in this work do not include observations of the velocity dispersion, 
the parameters <tro and R a are superfluous and may be ignored in fitting the galaxy. Likewise, 
only the major-axis surface brightness profile is used and therefore the disk scale-height 
parameter, Zd, may be fixed to a reasonable value. Finally, may be set to any value 
greater than 15kpc (i.e., beyond the outermost point of the observed rotation curve) and 
5rh may be ignored. The DF is thus specified by ten free parameters. 

Our set of model parameters must include the mass-to-light ratios of the disk and 
bulge. In general, the rotation curve fit for LSBs is improved by choosing a very large 
mass-to-light ratio for the disk, that is, by devising a galactic model that is disk-dominated 
in the inner regions. However, the required mass-to-light ratios are typically unphysical. 
Indeed, one can constrain mas s -to-light ratios u sing population synthesis models and galaxy 
colours (IBell fc de Jong! l200ll ; iBell et al.l 120031 ). In a Bayesian analysis such as MCMC, 
these constraints are implemented through prior probabilities for the mass-to-light ratios. 
For simplicity, we assume that these pri or probabilities follow a log-norm a l dist ribution. 
Using the B — R and B — V profiles from de Blok. van der Hulst. Sz Bothun Jl995h and the 
color-M/L relations from Bell et al. ( 20031 ) we find 



log((M/L) 



disk/ 1 



0.03 ±0.25 



(12) 
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and 



log Um/l) 



bulge 



0.2 ±0.25 



(13) 



The errors, which translate directly into the width of the prior probability distribution, are 
meant to incorporate uncertainties in the relations from iBell et al.l (120031 ). uncertainties in 
the colors, and differences in the M/L- values obtained by using either B — R or B — V 
colours. 

Two MCMC runs are conducted, one in which 7 is a free parameter and one in which 
7 is fixed to the NFW-value (i.e., 7 = 1). The surface brightness profile and rotation curve 
fits for a typical model from the first run a re shown in Figure HI Also shown is the fit 
found by iMcGaugh. Rubin, fc de Biota (120011 ) which assumes an exponential disk and does 
not include a bulge. Evidently, an excellent fit to the full surface brightness profile can be 
obtained provided both disk truncation and a bulge are included in the model. In Figure 
we show the probability distribution function for 7. Clearly, the data favor values of 7 
between and 0.8. 



The surface brightness profile and circular speed curve for a typical model from our 
MCMC analysis with 7 = 1 is shown in Figure El We find that the model r otation curve 



rises too quickly as compared with the data, in agreement with previous studies (IMoore 



Flores & Primack 


1994; 


McGaueh & de Blok 


1998 



1994 



4. MODELLING F-583 WITH A TRIAXIAL HALO 

Using the method outlined in Section 2, we transform the halo in one of our axisymmetric 
7 = 1 models. We begin by generating an N-body representation of the model with 500K 
particles for the halo, 400-fT particles for the disk, and 100K particles for the bulge. We 
produce two examples of models with triaxial halos: Model I where fl\ = 0.24 Gyr -1 and 
fa = 0.08 Gyr -1 (i.e., same choise of parameters as in our isolated halo example) and Model 
II where f3\ = 0.16 Gyr -1 and = — O.OSGyr^ 1 . Note that the artificial acceleration is 
applied only to the halo particles. 

First, consider Model I. Recall that in the example from Section 2 the short axis of the 
halo is aligned with the z-direction while the long axis is aligned with the x-direction. In a 
disk-bulge-halo system, the same choice of parameters leads to a rather mild deviation from 
axial symmetry since the intermediate axis is in the disk plane. 

For an isolated halo, the choice of parameters used in Model II amounts to a trivial 
interchange of the y and z coordinates. In the presence of the disk, the choice yields a halo 
model in which both the long and short axes are in the disk plane and hence the departure 
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from axial symmetry is very strong. 



In Figure we show contour plots of the disk surface density for Models I and II. As 
expected, departures from circular symmetry are more evident in Model I I Also, as expected, 
the lo ng axis of the disk is perpendicular to the long axis of the halo lHayashi fc Navarro 

(bood ). 



In Figure M we show the evolution of the halo axes ratios as a function of time for 
Models I and II. The evolution of the axis ratios in Model I is very similar to that found 
for the isolated halo in Section 2. By contrast, the influence of the disk is readily evident 
in Model II; the halo is flattened along the z-direction and somewhat rounder in the x — y 
plane than it would be in the absence of the disk. The end result is a halo that is prolate 
with approximate axial symmetry about the x-axis. 

In Figure [9] we show the rotation curves for Models I and II as calculated along a slit 
placed on the major axis of the disk. In both experiments, the maximum rotation speeds 
decreases relative to their initial values by about lOkms -1 . One might imagine an iterative 
procedure in which, given these results, one adjusts the initial model so that the final system 
better reproduces the data. 

Next, we consider the change in shape of the rotation curve produced by the defor- 
mation of the halo. In both experiments, the rotation curve rises more slowly than in the 
initial, axisymmetric model. Figure M shows the logarithmic slope of the rotation curve and 
illustrates this point quantitatively. We can also quantify the change in shape of the rotation 
curve by considering the fitting formula 



v(r) = v - 



1 



1 + x c 



A/a 



(14) 



where x 



Courteau 


(1997 


). 


( 


Courteau! ( 


1997) 



for our purposes, this form will suffice (see, for example, (lHayashi et al.l 12004 ) .)) A cored- 
isothermal sphere yi elds a rotation c urve with a > 2 wh ile an NFW-halo yields a rotation 
curve with a ~ 0.6 (ICourteaulll997l ; lHayashi et al.ll2004l ). We find the following values for 



a: observed rotation curve - 5.35; initial, axisymmetric model - 0.91; Model I - 1.94; and 
Model II - 3.57. Clearly, Model II, where departures from axial symmetry are strongest, 
comes closest to reproducing the shape of the rotation curve. 
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SUMMARY AND CONCLUSION 



The adiabatic squeezing method lHolley-Bockelmann et al.l (120011 ) produces triaxial halos 
that have shrunk in size and therefore requires that the positions and velocities of the particles 
be rescaled. This awkward step precludes the technique from being applied to compound 
systems. Our approach avoids this problem by using an energy-conserving artificial force to 
deform the halos. 

Our analysis of the LSB galaxy F568-3 begins with a discussion of axisymmetric models. 
We attempt to fit both photometric and kinematic observations using Bayesian statistics and 
the MCMC method. Our excellent fit of the surface brightness profile requires a bulge and 
disk truncation, neither of which were included in previous studies. As for the rotation curve, 
we find that constant density cores do better than density cusps in agreement with earlier 
studies of LSBs. 

The second stage of our analysis is to deform the halo of a compound system. In 
agreement with lHayashi fc Navarro! (120061 ). we show that the rotation curve of F568-3 may 
indicate the presence of a triaxial halo rather than a pr oblem with the standard CDM mode l 
of structure formation. lHayashi fc Navarro! (120061 ) and lHayashi. Navarro, fc Springell (120071 ) 
construct rotation curves by finding closed orbits in the gravitational potential of a triaxial 
halo. We calculate the rotation curves by making pseudo-observations of a disk that is 
self-consistently embedded in a dark halo. 

There are two improvements that will add a further level of realism to the analysis: the 
inclusion of a gas disk in the galactic models and an iterative scheme whereby the initial 
model and artificial acceleration parameters are adjusted so that the final model fits the data 
in detail. These improvements will be considered in a future publication. 

Our triaxial models have a wide range of applications. For example, they can be used 
to study the effect a non-spherical halo has on the morphology of tidal streams from satellite 
galaxies and flaring and warping of the gas disk. The method can also be applied to bulges 
where departures from axial symmetry are thought to be important. 



It is a pleasure to thank J. Bailin, S. Courteau, J. Dubinski, S. McGaugh, and D. 
Puglielli for useful conversations. This work was supported by a grant from the Natural 
Sciences and Engineering Research Council of Canada. 
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Fig. 1. — Surface density contours along the three principle axes for the triaxial isolated halo 
constructed in Section 2. Spacing between solid contours is 1 dex. 
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Fig. 2. — Axes ratios as a function of time. Solid curves show b; dashed curves show c. Red, 
blue, and green curves are for, respectively, the inner one third of the particles, the inner 
two thirds of the particles, and the entire halo. 
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Fig. 3. — Differential mass profile, dM/dr oc r 2 p as a function of spherical radius r. Line 
types are: analytic profile - solid black curve; initial profile derived from the N-body dis- 
tribution - dot-dashed magenta curve; profile at 6 Gyr - dotted red curve; profile at 12 Gyr 
- dashed blue curve. The profile for the control experiment (no artificial acceleration) at 
12 Gyr is shown by the long-dashed green curve. The straight solid black line corresponds to 
r 2 p oc r or p(r) oc r _1 . Lower panel gives the fractional difference between the 4 measured 
profiles and the the analytic expression (i.e., (p meaS ured - Pcxact) //Wet- 



Table 1. Invariances of Equations 1-3 



P Coordinates Axis ratios 



Pi — > Pi Pi^> Pi x y y — > rr 2-^2 6— > 1/6 c — > c/6 
Pi ^ Pi — Pi Pi — > — /?2 x ^ :r y ^ z z — > rr 6 — > c c — > 6 
Pi^ P2- Pi P2 — >■ -A 2/^2 6 — > 1/c c^b/c 



Table 2. Axis Ratios for Various Choices of /3i and /3 2 



Isolated Halo Composite Model 



Model 


A (Gyr- 1 ) 


A (Gyr- 1 ) 


b 


c 


b 


c 


la 


0.12 


0.04 


0.93 


0.86 


0.95 


0.85 


lb 


0.08 


-0.04 


0.86 


0.93 


0.88 


0.90 


lc 


-0.08 


-0.12 


0.93 


1.07 


0.93 


1.06 


2a 


0.16 


0.16 


1.00 


0.78 


1.01 


0.77 


2b 


0.0 


-0.16 


0.78 


1.00 


0.80 


0.99 


3a 


0.24 


0.08 


0.89 


0.73 


0.92 


0.73 


3b 


0.16 


-0.08 


0.73 


0.89 


0.75 


0.87 


3c 


-0.16 


-0.24 


0.82 


1.12 


0.82 


1.08 
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Fig. 4. — Comparison of data with predictions for a typical model from the MCMC run 
wher e 7 is a free parameter. Top panel sho ws surface brightness profile. Observations 
from Ide Blok. van der Hulst. fc Bothunl (119951 ) are indicated by magneta dots. Line types 
are as follows: total model surface brightness profile - solid black curve; disk contribution 
- blue dotted curve; bulge contribution - red dashed curve; exponential disk model from 



de Blok. van der Hulst. fc Bothunl (119951) - yellow long-dashed curve. Botto m panel shows 



the rotation curve. Observations from iMcGaugh. Rubin, fc de Blok! (120011 ) are indicated 
by magneta dots. Total model - solid black curve; disk contribution - blue dotted curve; 
bulge contribution - red dashed curve; halo contribution - green long-dashed curve; gas 
contribution - yellow dot-dashed curve. 
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Fig. 5. — Probability distribution function for 7 from the first MCMC where 7 is a free 
parameter. 
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Fig. 6. — Surface brightness profile and circular speed curve for a typical mo 
chain where 7 = 1. Line types are the same as in Figure HI 



del from MCMC 
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Fig. 7. — Superposition of surface density contours and projected particle distribution for 
Models A (top) and B (bottom). 
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Fig. 8. — Halo axis ratios as a function of time for Models A (top) and B (bottom). Line 
types are the same as in Figures [2j 
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Fig. 9. — Observed and model rotation curves. Top pa nel shows the rotation curves; bottom 
panel shows the logarithmic slope. Observations from iMcGaugh. Rubin, fc de Blokl (120011 ) 
are shown as red dots. Line types are as follows: black dotted line - rotation curve for the 
initial model as derived from the potential; blue dashed line - rotation for the initial model 
as derived by measuring bulk motion of the disk stars; green long-dashed line - rotation 
curve for Model I; magenta solid line - rotation curve for Model II. 



